MODULE smm
  USE GLOBVAR
  USE NRUTIL
  USE PROBABILITY
  USE MATRIX
  USE STATISTICS
  USE random
  IMPLICIT NONE
  
    CONTAINS
   REAL(8) FUNCTION mom(thetain) 
   !SUBROUTINE mom(theta_ini,beta_f,check,N_int,sc)
        IMPLICIT NONE
        REAL(8), INTENT(IN):: thetain(:) !initial guess

        REAL(8):: Ls(N_inds,N_testl),Ep(N_inds,N_testl),as,pen,pen2
        REAL(8):: L0s(N_inds),eta(N_inds,N_testl),e_L0(n_inds), e_eta(n_inds),mom_simu(N_mom),dif(N_mom,1),inters(n_inds,n_levell,20,3),interls(n_inds,n_levell,3)!, latent(N_inds*N_testl,N_test),ps(N_testl,N_test)
        REAL(8),ALLOCATABLE:: AVSE(:),cors(:,:)
        ! parameters
        REAL(8):: beta_mu(6,1),beta_eta(n_x_eta,1),var_l,Lbar(N_testl),var_e(N_testl), etav,delta(N_levell),beta_level(n_l,1),beta_y(2,1)
        INTEGER:: Ds(N_inds,N_testl),mks(N_testl,N_test),nks(N_testl,N_test),indexl(n_testl,2), artls(n_inds,11,20),duls(n_inds,n_levell),dus(n_levell),inds(N_inds,n_levell) 
        INTEGER:: i,j,k,jj,kk,m,t,n,l,mmm
        
        !!!!!!!!!!!!!!!!!!Simulate Moments
        CALL SET_SEED(1438914432,727629,19890414,111111)
        beta_mu(1:n_x_l0,1)=thetain(1:n_x_l0)                                              ! Initial Latent skill
        DO i=2,4
          beta_mu(i,1)=1.40d0**(thetain(i))/(10.0d0*(1.0d0+1.40d0**(thetain(i))))
        END DO
        beta_mu(5,1)=1.0d0-1.40d0**(thetain(5))
        
        var_l=4.0d0*(1.4d0**thetain(n_x_l0+1)/(1.0d0+1.4d0**thetain(n_x_l0+1)))     ! variance of the shoch of initial latent skill
                Lbar=1.0d0
        Lbar(1)=0.0d0
        Lbar(3:N_levell) =thetain(n_x_l0+2:n_x_l0+N_levell-1)                ! mimimum requriement for each difficulty level
        var_e=1.0d0
        DO i=3,n_levell
          var_e(i)=1.4d0**thetain(n_x_l0+n_levell-1+i-2)     ! variance of the shock at each difficulty level
        END DO
        beta_eta(1:n_x_eta,1)=thetain(n_x_l0+2*N_levell-3+1:n_x_l0+2*N_levell+n_x_eta-3)           ! x for eta
        etav=1.40d0**thetain(n_x_l0+2*N_levell+n_x_eta-2)                           ! variance of shock for eta
        delta=1.0d0
        DO i=3,N_levell
            delta(i)=6*(1.10d0**thetain(n_x_l0+2*N_levell+n_x_eta-2+i-2))/(1.0d0+1.10d0**thetain(n_x_l0+2*N_levell+n_x_eta-2+i-2))  ! delta for each level
        END DO
        beta_level=0.0d0
        beta_level(n_levell:n_l,1)=thetain(n_x_l0+2*N_levell+n_x_eta-2+n_levell-2+1:n_x_l0+2*N_levell+n_x_eta-2+n_levell-2+n_levell-1)
                beta_y=0.0d0
                     beta_y(1:2,1)=thetain(n_x_l0+2*N_levell+n_x_eta-2+n_levell-2+n_levell-1+1:n_x_l0+2*N_levell+n_x_eta-2+n_levell-2+n_levell-1+2)

        ! generate shock terms
        
        DO i=1,N_inds
            e_L0(i) =Sample_NORMAL(0.0d0,var_l)
            k=i-INT((i-1)/N_IND)*N_ind
            e_eta(i)=Sample_Normal(intermean(k),etav)
        END DO
        
        ! GENERATE INITIAL LATENT VALUES AND 
        DO i=1,N_inds
            k=i-INT((i-1)/N_IND)*N_ind
            L0s(i)=1.40d0**(DOT_PRODUCT(L0(k,:),beta_mu(:,1))+e_L0(i))
        END DO
        
       ! SIMULATE LATENT SKILL PROCESS ! HERE ONLY FOR LANGUAGE SKILLS 
        
        DO i=1,N_inds
            k=i-INT((i-1)/N_IND)*N_ind
            DO j=1,N_testl
                IF (j<firstl(k,1)) THEN
                    Ls(i,j)=-99.0d0
                ELSE IF (j==firstl(k,1) .and. m_l(k,j)==1) THEN
                    IF (x_testl(k,j,1)<-90) THEN
                        eta(i,j)=1.40d0**(e_eta(i))/((1.0d0+1.40d0**(e_eta(i))))
                    ELSE
                        eta(i,j)=1.40d0**(DOT_PRODUCT(x_testl(k,j,:),beta_eta(:,1))+e_eta(i))/((1.0d0+1.40d0**(DOT_PRODUCT(x_testl(k,j,:),beta_eta(:,1))+e_eta(i))))
                    END IF
                    
                    kk=levell(j)
68                    Ls(i,j)=L0s(i)*(exp(delta(kk)*eta(i,j)+Sample_Normal(0.0d0,var_e(kk))+beta_y(1,1)*monthl(j)+beta_y(2,1)*DBLE(timel(j))))
                    !IF (Ls(i,j)<0.0d0) GO TO 68 
                ELSE IF (j==firstl(k,1) .and. m_l(k,j)==-1) THEN
                    Ls(i,j)=-99.0d0
                ELSE IF (j>firstl(k,1) .and. Ls(i,j-1)<0.0d0 .and. m_l(k,j)==-1) THEN
                    Ls(i,j)=-99.0d0
                ELSE IF (j>firstl(k,1) .and. Ls(i,j-1)<0.0d0 .and. m_l(k,j)==1) THEN
                    IF (x_testl(k,j,1)<-90) THEN
                        eta(i,j)=1.40d0**(e_eta(i))/((1.0d0+1.40d0**(e_eta(i))))
                    ELSE
                        eta(i,j)=1.40d0**(DOT_PRODUCT(x_testl(k,j,:),beta_eta(:,1))+e_eta(i))/((1.0d0+1.40d0**(DOT_PRODUCT(x_testl(k,j,:),beta_eta(:,1))+e_eta(i))))
                    END IF
                    kk=levell(j)
77                    Ls(i,j)=L0s(i)*(exp(delta(kk)*eta(i,j)+Sample_Normal(0.0d0,var_e(kk))+beta_y(1,1)*monthl(j)+beta_y(2,1)*DBLE(timel(j))))
                    !IF (Ls(i,j)<0.0d0) GO TO 77
                ELSE IF (j>firstl(k,1) .and. Ls(i,j-1)>=0.0d0 .and. m_l(k,j)==-1 .and. (levell(j)==levell(j-1))) THEN
                    Ls(i,j)=Ls(i,j-1)
                ELSE IF (j>firstl(k,1) .and. Ls(i,j-1)>=0.0d0 .and. m_l(k,j)==-1 .and. (levell(j) .ne. levell(j-1))) THEN
                    kk=levell(j)
                    IF (kk>0) THEN
                       Ls(i,j)=(beta_level(kk-1,1)*1.0d0+beta_level(n_levell-1+kk-1,1)*Ls(i,j-1))  
                    ELSE 
                        Ls(i,j)=Ls(i,j-1)
                    END IF
                ELSE IF (j>firstl(k,1) .and. Ls(i,j-1)>=0.0d0 .and. m_l(k,j)==1 .and. (levell(j)==levell(j-1))) THEN
                    IF (x_testl(k,j,1)<-90) THEN
                        eta(i,j)=1.40d0**(e_eta(i))/((1.0d0+1.40d0**(e_eta(i))))
                    ELSE
                        eta(i,j)=1.40d0**(DOT_PRODUCT(x_testl(k,j,:),beta_eta(:,1))+e_eta(i))/((1.0d0+1.40d0**(DOT_PRODUCT(x_testl(k,j,:),beta_eta(:,1))+e_eta(i))))
                    END IF
                    kk=levell(j)
84                  Ls(i,j)=Ls(i,j-1)*(exp(delta(kk)*eta(i,j)+Sample_Normal(0.0d0,var_e(kk))+beta_y(1,1)*monthl(j)+beta_y(2,1)*DBLE(timel(j))))
                    !IF (Ls(i,j)<0.0d0) GO TO 84
                ELSE IF (j>firstl(k,1) .and. Ls(i,j-1)>=0.0d0 .and. m_l(k,j)==1 .and. (levell(j) .ne. levell(j-1))) THEN
                    IF (x_testl(k,j,1)<-90) THEN
                        eta(i,j)=1.40d0**(e_eta(i))/((1.0d0+1.40d0**(e_eta(i))))
                    ELSE
                        eta(i,j)=1.40d0**(DOT_PRODUCT(x_testl(k,j,:),beta_eta(:,1))+e_eta(i))/((1.0d0+1.40d0**(DOT_PRODUCT(x_testl(k,j,:),beta_eta(:,1))+e_eta(i))))
                    END IF
                    kk=levell(j)
92                  Ls(i,j)=((beta_level(kk-1,1)*1.0d0+beta_level(n_levell-1+kk-1,1)*Ls(i,j-1)))*(exp(delta(kk)*eta(i,j)+Sample_Normal(0.0d0,var_e(kk))+beta_y(1,1)*monthl(j)+beta_y(2,1)*DBLE(timel(j))))
                    !IF (Ls(i,j)<0.0d0) GO TO 92    
                END IF
            END DO
        END DO
        
        ! GENERATE TASK PERFORMANCE !
        
        DO i=1,N_inds
            k=i-INT((i-1)/N_IND)*N_ind
            DO j=1,N_testl
                kk=levell(j)
                IF (m_l(k,j)==-1) THEN
                    DS(i,j)=-99
                ELSE IF (Ls(i,j)>=SUM(Lbar(1:kk)) .and. m_l(k,j)==1) THEN
                    DS(i,j)=1
                ELSE IF (Ls(i,j)<SUM(Lbar(1:kk)) .and. m_l(k,j)==1) THEN
                    DS(i,j)=0
                END IF
            END DO
        END DO
                
         !!!!!!!!!!!!!!!!!From here !!!!!!!!!!!!!!!!!!!!!!!
        
    !!! language tasks!!!!
    indexl=0
    DO i=1,N_inds
      DO j=1,N_testl
         IF (ds(i,j)>-90) THEN
             indexl(j,1)=indexl(j,1)+1
             IF (ds(i,j)>0.5) indexl(j,2)=indexl(j,2)+1
         END IF         
      END DO
    END DO
    
    DO i=1,N_testl
        IF (indexl(i,1)>0.5) THEN
            mom_simu(i)=DBLE(indexl(i,2))/DBLE(indexl(i,1)) ! PASSING RATE FOR EACH LANGUAGE TASK
        ELSE 
            mom_simu(i)=-99
        END IF
    END DO
    
    !!!! generate artl variables !!!!
    artls=-1
    DO i=1,N_inds
        DO j=1,N_testl
            IF (ds(i,j)>-90) THEN
                kk=levell(j)
                k=i-INT((i-1)/N_IND)*N_ind
                m=repl(k,j)
                artls(i,kk,m)=ds(i,j)
                inters(i,kk,m,1:3)=x_testl(k,j,1:3)
            END IF
        END DO
    END DO
    
    !!!! passing rate at each level!!!!

    DO j=N_testl+1,N_testl+N_levell
       m=0
       n=0 
       DO i=1,N_inds
           DO k=1,20
               IF (artls(i,j-N_testl,k)>-0.5) THEN
                   m=m+1
                   IF (artls(i,j-N_testl,k)==1) n=n+1
               END IF
           END DO 
       END DO
       IF (m>0) THEN
           mom_simu(j)=DBLE(n)/DBLE(m)
       ELSE
           mom_simu(j)=-99
       END IF
       
    END DO
    
    !!!! correlation across levels !!!!
    
    DO j=2,N_levell-1
      DO t=1,5
        DO k=1,5
          m=0
          n=0 
          DO i=1,N_inds 
             IF (artls(i,j,t)==1 .and. artls(i,j+1,k)>-0.5) THEN    ! the probability of passing the kth task at level j+1 conditional on passing the tth task at level j
                 m=m+1
                 IF (artls(i,j,t)==1 .and. artls(i,j+1,k)==1) n=n+1
             END IF
          END DO
          
          IF (m>0) THEN
             mom_simu(N_testl+N_levell+(j-2)*25+(t-1)*5+k)=DBLE(n)/DBLE(m)
          ELSE
             mom_simu(N_testl+N_levell+(j-2)*25+(t-1)*5+k)=-99
          END IF
        END DO
      END DO
   END DO        
        
  DO j=2,N_levell
        l=0
      DO t=1,4
        DO k=t+1,5
          m=0
          n=0 
          l=l+1
          DO i=1,N_inds 
             IF (artls(i,j,t)==1 .and. artls(i,j,k)>-0.5) THEN    ! the probability of passing the kth task at level j+1 conditional on passing the tth task at level j
                 m=m+1
                 IF (artls(i,j,t)==1 .and. artls(i,j,k)==1) n=n+1
             END IF
          END DO
          
          IF (m>0) THEN
             mom_simu(25*(N_levell-2)+N_testl+N_levell+(j-2)*10+l)=DBLE(n)/DBLE(m)
          ELSE
             mom_simu(25*(N_levell-2)+N_testl+N_levell+(j-2)*10+l)=-99
          END IF
        END DO
      END DO
  END DO       
  
  DO j=2,N_levell-2
      DO t=1,5
        DO k=1,5
          m=0
          n=0 
          DO i=1,N_inds 
             IF (artls(i,j,t)==1 .and. artls(i,j+2,k)>-0.5) THEN    ! the probability of passing the kth task at level j+1 conditional on passing the tth task at level j
                 m=m+1
                 IF (artls(i,j,t)==1 .and. artls(i,j+2,k)==1) n=n+1
             END IF
          END DO
          
          IF (m>0) THEN
             mom_simu(25*(N_levell-2)+N_testl+N_levell+(N_levell-1)*10+(j-2)*25+(t-1)*5+k)=DBLE(n)/DBLE(m)
          ELSE
             mom_simu(25*(N_levell-2)+N_testl+N_levell+(N_levell-1)*10+(j-2)*25+(t-1)*5+k)=-99
          END IF
        END DO
      END DO
  END DO
  
  indexl=0   ! task passing rate for new enrolled children
    DO i=1,N_inds
        k=i-INT((i-1)/N_IND)*N_ind
      DO j=1,N_testl
          as=DBLE(firstl(k,2)-monthl(j)-DBLE(timel(j))/4.0d0)
         IF (ds(i,j)>-90 .AND. as>-1.25) THEN
             indexl(j,1)=indexl(j,1)+1
             IF (ds(i,j)>0.5) indexl(j,2)=indexl(j,2)+1
         END IF         
      END DO
    END DO
    
    DO i=1,N_testl
        IF (indexl(i,1)>0.5) THEN
            mom_simu(25*(2*N_levell-5)+N_testl+N_levell+(N_levell-1)*10+i)=DBLE(indexl(i,2))/DBLE(indexl(i,1)) ! PASSING RATE FOR EACH LANGUAGE TASK
        ELSE 
            mom_simu(25*(2*N_levell-5)+N_testl+N_levell+(N_levell-1)*10+i)=-99
        END IF
    END DO
    
    
    indexl=0   ! task passing rate for new enrolled children
    DO i=1,N_inds
        k=i-INT((i-1)/N_IND)*N_ind
      DO j=1,N_testl
          as=DBLE(firstl(k,2)-monthl(j)-DBLE(timel(j))/4.0d0)
         IF (ds(i,j)>-90 .AND. as<-1.25) THEN
             indexl(j,1)=indexl(j,1)+1
             IF (ds(i,j)>0.5) indexl(j,2)=indexl(j,2)+1
         END IF         
      END DO
    END DO
    
    DO i=1,N_testl
        IF (indexl(i,1)>0.5) THEN
            mom_simu(25*(2*N_levell-5)+2*N_testl+N_levell+(N_levell-1)*10+i)=DBLE(indexl(i,2))/DBLE(indexl(i,1)) ! PASSING RATE FOR EACH LANGUAGE TASK
        ELSE 
            mom_simu(25*(2*N_levell-5)+2*N_testl+N_levell+(N_levell-1)*10+i)=-99
        END IF
    END DO
    
       l=0
    DO j=2,N_levell
      DO t=1,5
          m=0
          n=0
          l=l+1
          DO i=1,N_inds 
             IF (artls(i,j,t)>-0.5) THEN   
                 m=m+1
                 IF (artls(i,j,t)==1) n=n+1
             END IF
          END DO
          
          IF (m>0) THEN
             mom_simu(25*(N_levell-2)+3*N_testl+N_levell+(N_levell-1)*10+25*(N_levell-3)+l)=DBLE(n)/DBLE(m)
          ELSE
             mom_simu(25*(N_levell-2)+3*N_testl+N_levell+(N_levell-1)*10+25*(N_levell-3)+l)=-99
          END IF
      END DO
    END DO
    
    mmm=25*(N_levell-2)+3*N_testl+N_levell+(N_levell-1)*10+25*(N_levell-3)+5*(N_levell-1)
    
    DULs=-99
    inds=-99
    DO i=1,N_inds
        DO k=1,N_levell
            DO j=1,20
                IF (artls(i,k,j)==1 .and. j==1) THEN
                    duls(i,k)=1
                    inds(i,k)=1
                ELSE IF (artls(i,k,j)==0 .and. j==1) THEN
                    duls(i,k)=j
                    inds(i,k)=0
                ELSE IF (j>1 ) THEN
                    IF( artls(i,k,j-1)==0 .and. artls(i,k,j)==0 .and. inds(i,k)==0) THEN
                       duls(i,k)=j
                       inds(i,k)=0
                       IF (j<20 .and. artls(i,k,j+1)==-1) THEN
                           duls(i,k)=j+1
                           inds(i,k)=1
                       END IF  
                    ELSE IF (j>1 .and. artls(i,k,j-1)==0 .and. artls(i,k,j)==0 .and. inds(i,k)==1) THEN
                        duls(i,k)=duls(i,k)
                        inds(i,k)=1
                    ELSE IF (j>1 .and. artls(i,k,j-1)==0 .and. artls(i,k,j)==1 .and. inds(i,k)==0) THEN
                        duls(i,k)=j
                        inds(i,k)=1
                    ELSE IF (j>1 .and. artls(i,k,j-1)==0 .and. artls(i,k,j)==1 .and. inds(i,k)==1) THEN
                        duls(i,k)=duls(i,k)
                        inds(i,k)=1   
                    END IF
                END IF
                
            END DO           
                
        END DO
    END DO
    
    DO k=1,n_levell
        m=0
        n=0
        DO i=1,n_inds
            IF (duls(i,k)>-90) THEN
               m=m+1
               n=n+duls(i,k)
            END IF
        END DO
        IF (m>0) THEN
          mom_simu(mmm+k)=DBLE(n)/DBLE(m)
        ELSE 
           mom_simu(mmm+k)=-99
        END IF
    END DO
    
    mmm=25*(N_levell-2)+3*N_testl+N_levell+(N_levell-1)*10+25*(N_levell-3)+5*(N_levell-1)+n_levell
    
        dus=0
    DO k=1,n_levell
        DO i=1,N_inds
            IF (duls(i,k)>-90) THEN
                dus(k)=dus(k)+1
            END IF
        END DO
    END DO
    DO k=1,n_levell
        allocate(avse(dus(k)))
        m=0
        DO i=1,n_inds
            IF (duls(i,k)>-90) THEN
                m=m+1
                avse(m)=duls(i,k)
            END IF
        END DO
        mom_simu(mmm+k)=STDEV_V(AVSE)
        DEALLOCATE(avse)
    END DO
     mmm=25*(N_levell-2)+3*N_testl+N_levell+(N_levell-1)*10+25*(N_levell-3)+5*(N_levell-1)+2*n_levell

    INDs=0
    interls=0.0d0
    DO k=2,N_levell
        DO i=1,N_inds
          DO j=1,20
              IF (inters(i,k,j,1)>-90) THEN
                  inds(i,k)=inds(i,k)+1
                  interls(i,k,1)=interls(i,k,1)+inters(i,k,j,1)
                  interls(i,k,2)=interls(i,k,2)+inters(i,k,j,2)
                  interls(i,k,3)=interls(i,k,3)+inters(i,k,j,3)
              END IF
          END DO
          IF (inds(i,k)>0) THEN
               interls(i,k,1)=DBLE(interls(i,k,1))/DBLE(inds(i,k))
               interls(i,k,2)=DBLE(interls(i,k,2))/DBLE(inds(i,k))
               interls(i,k,3)=DBLE(interls(i,k,3))/DBLE(inds(i,k))
          END IF
        END DO
    END DO
    
    dus=0
    DO k=1,n_levell
        DO i=1,N_inds
            IF (inds(i,k)>0 .AND. duls(i,k)>-90) THEN
            dus(k)=dus(k)+1          
            END IF
        END DO
    END DO
    
    DO k=1,n_levell
        allocate(cors(dus(k),4))
        m=0
        DO i=1,n_inds
             IF (inds(i,k)>0 .AND. duls(i,k)>-90) THEN
                m=m+1
                cors(m,1)=interls(i,k,1)
                cors(m,2)=interls(i,k,2)
                cors(m,3)=interls(i,k,3)
                cors(m,4)=duls(i,k)
             END IF
        END DO
        IF (k>1) THEN
        mom_simu(mmm+3*(k-1)+1)=Correlation(cors(:,1),cors(:,4))
        mom_simu(mmm+3*(k-1)+2)=Correlation(cors(:,2),cors(:,4))
        mom_simu(mmm+3*(k-1)+3)=Correlation(cors(:,3),cors(:,4))
        ELSE 
            mom_simu(mmm+3*(k-1)+1)=-99
            mom_simu(mmm+3*(k-1)+2)=-99
            mom_simu(mmm+3*(k-1)+3)=-99
        END IF
        
            
        DEALLOCATE(cors)
    END DO
    mmm=25*(N_levell-2)+3*N_testl+N_levell+(N_levell-1)*10+25*(N_levell-3)+5*(N_levell-1)+5*n_levell
    
    
    
    !    
    !     !!!!!!!!!!!!! Now finish simulate the fake data and start to generat the moments !!!!!!!!!!!!!!!!!!!!!!!!!
    !     
    
    DO i=1,N_mom
        IF (ABS(moment(i))<1.0d-10) THEN
            dif(i,1)=((mom_simu(i)-moment(i)))**2.0D0
        ELSE 

                        dif(i,1)=((DBLE((mom_simu(i)-moment(i))))**2.0D0)
        END IF 
        IF (i<=N_testl+N_levell) THEN
           dif(i,1)=1.0d0*dif(i,1)
        END IF 

    END DO
    

    mom=SUM(dif)
    
    
    
                OPEN(6424,file="dif.out")
                  DO i = 1,N_mom
                    WRITE(6424,'(200F16.8)') SQRT(dif(i,1)),mom_simu(i),moment(i)
                  END DO
                CLOSE(6424)
            
   
   END FUNCTION mom
    
       
    
END MODULE smm 